Hierarchical theory of quantum dissipation: Partial fraction decomposition scheme 
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We propose a partial fraction decomposition scheme to the construction of hierarchical equations 
of motion theory for bosonic quantum dissipation systems. The expansion of Bose-Einstein function 
in this scheme shows similar properties as it applies for Fermi function. The performance of the 
resulting quantum dissipation theory is exemplified with spin-boson systems. In all cases we have 
tested the new theory performs much better, about an order of magnitude faster, than the best 
available conventional theory based on Matsubara spectral decomposition scheme. 



I. INTRODUCTION 

It is well established now that path integral influence 
functional formalism [3, 0, El of quantum dissipation the- 
ory (QDT) can be reformulated in terms of hierarchical 
equations of motion (HEOM) 0, 1, i, 0, i, B 0, HH, O • 
Compared to its dynamics equivalent path integral influ- 
ence functional formalism, the HEOM has the advantage 
in both numerical efficiency and applications to various 
systems. Moreover, the initial system-bath coupling that 
is not contained in the original path integral formalism 
can now be accounted for via proper initial conditions to 
HEOM [H[I3]. As an exact formalism, HEOM is also 
nonperturbative. It treats the combined effects of sys- 
tem anharmonicity, system-bath coupling strength, and 
multiple memory time scales on the reduced system dy- 
namics. 

The specific HEOM-QDT construction depends how- 
ever on the way of decomposing bath correlation func- 
tion into its memory components. Different decompo- 
sition schemes are mathematically equivalent, but have 
different numerical performance. To illustrate this is- 
sue, we consider only the single dissipation mode case, 
in which the system-bath coupling Hamiltonian assumes 
Hsb = —QF B . The system operator Q here is also 
called the dissipation mode, through which the gener- 
alized Langevin force F B (t) = e i- h Bt/hp^ e -ih B t/h ac ^ s Qn 

the system. As exact theory is concerned, the bath cor- 
relation function C(t — r) = (F B (t)F B (T)) is related to 
the bath spectral density function J(tv) via 



C(t) 



duj e 



1 - ' 



(1) 



This is the fluctuation-dissipation theorem of bosonic 
canonical ensembles, where (3 = h/(k,BT) denotes the 
inverse temperature. 
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The common strategy of memory decomposition used 
for HEOM construction is based on the Matsubara series 
expansion of Bose (also called Bose-Einstein) function in 
Eq. |T]). It assumes (setting x = (3lo) 



l-e~ x ~ 2 + x + ^ 



m—1 



l 



l 



x + i2i:m x — ilixm 



(2) 



This is the Matsubara spectral decomposition (MSD) 
scheme. It is exact when N — > oo. It resolves the memory 
contents of bath correlation function with Matsubara fre- 
quencies (27rm//3; m > 1), together with the poles of J{z) 
in lower plane (Imz < 0). In practical use, the residue 
6C(t), due to the deviation of Eq. @ from the exact one, 
is treated by white noise [13(. The resulting HEOM for- 
malism has been applied to electron/population transfer 
and optical spectroscopy systems [H, HH, [13] • 

In this work, we explore the partial fraction decompo- 
sition (PFD) scheme. This scheme was originally pro- 
posed for Fermi function in the study of electronic dy- 
namics/structure systems (l8|. [11 |20j [2lj [22] [23|. |24|. 
The motivation behind is that MSD is well-known to 
be of slow convergence. The PFD of Bose function will 
be carried out in a similar manner as that proposed re- 
cently by Croy and Saalmann for the Fermi function 
[2~il ]. We present the PFD results on Bose function in 
Sec. |H]and the corresponding HEOM-QDT construction 
in Sec. IIIII Their derivations are given in Appendix lAl 
and Appendix [Bl respectively. 

Numerical performance of the new HEOM-QDT will 
be exemplified with a model spin-boson Drude dissi- 
pation system in Sec. HVl Included are also compar- 
isons with the best available MSD-based conventional 
HEOM theory. The latter is implemented with the well- 
established Markovian residue correction method @, 0] 
that has improved significantly the performance of MSD- 
based HEOM. Nevertheless, the PFD-based HEOM the- 
ory without residue correction still performs much better, 
about 20 times faster, than the best available MSD-based 
approach. Finally we conclude this work. 
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FIG. 1: Poles of the PFD expansion (solid circle) at the ex- 
pansion order N=2, 4, 8, 16. Included for comparison are also 
the Matsubara poles (cross). Poles are complex-conjugate 
paired, and only those in lower plane are shown. 




FIG. 2: The difference <5/jv(x) = /be (a?) — /n(x), between 
the exact Bose-Einstein function and its PFD expansion ap- 
proximation [Eq. 0], as the function of x at 7V=2 (black), 4 
(red), 8 (blue), and 16 (orange). The inset is for <S/jj(a;) as 
function of TV, at x=5 (black) and 50 (red). The Matsubara 
counterparts are also shown by dashed curves. 



II. PARTIAL FRACTION DECOMPOSITION 
SCHEME 

A. Bose— Einstein function in decomposition 

The PFD scheme starts with the identity 

1 11 cosh(:c/2) 

l-e~ x 2 2 smh{x/2) ' [ > 

followed by Taylor's expansions of the numerator func- 
tion cosh(a;/2) to the (27V) th -order and the denominator 
sinh(x/2) to the (27V + l) th -orders, respectively. PFD 
leads to Bose function the expression (cf. Appendix [A"|l 

1 _ 1 1 y> / 1 1 \ 

1-e" 35 ~ 2 + x + ^U + 2 v /^ + x-2y/^)' 

The involving parameters {£.,■; j = 1, • • • TV} are the roots 

of the denominator polynomial, J2n=o^ n /(^ n + -0' = 0- 
They can be determined as the eigenvalues of an TV x TV 
matrix whose elements are 

A mn = 2m(2m + l)S m+hn - 2N{2N + l)S mN . (5) 

Note that the Fermi function counterpart is similar [24| , 
just replacing +1 with —1 in both parentheses in Eq. ([5]). 
The derivation of the above formalism is detailed in Ap- 
pendix El 

Figure [TJ depicts the poles Xj = ±2y^J of the PFD 
scheme at different orders. Shown are only those N 
poles in lower plane, and their complex conjugates (in 
upper plane) are implied. The distribution pattern of 
the total 27V poles is similar to that of Fermi function 
[23 | . There are not only pure imaginary poles, which are 
mostly close to the Matsubara poles, but also complex 



poles with nonzero real parts. This is right the feature 
for the efficiency of PFD scheme. 

Figure [2] depicts the deviation of the approximations 
from the exact Bose function. The curves there can be 
used to estimate the required order of TV for the PFD- 
HEOM dynamics, by considering the effective system and 
bath frequencies. It is clearly seen that the PFD expan- 
sion well overlaps with the exact result, within the range 
\x\ < 47V (i.e., \oj\ < ANk B T/K), as estimated, although 
for the case of Fermi function by Croy and Saalmann 
recently [2"3 |. 



B. Correlation function in decomposition 

For the later construction of HEOM formalism, we 
shall expand the bath correlation function C(t) in an 
exponential series @. This can be achieved via contour 
integration method. It recasts Eq. ([TJ) via analytic con- 
tinuation as 

C{t)= 1 -jdze-^ T ^ Fz , (6) 

and evaluates it by Cauchy's residue theorem. The con- 
tour of integration encloses the lower-half plane for the 
required C(t > 0). Denote 

C(t) = C (t)+C B (t), (7) 

for the contributions from the poles of spectral den- 
sity J(z) and Bose function, respectively. Note that 
the antisymmetry property of bosonic spectral density, 
J(— u>) — —J(ui), reads now 

J{-z*) = -J*(z). (8) 
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Consider first the Bose function contribution C-a(t). 
The total N PFD poles in lower plane are either pure 
imaginary or true complex numbers in pairs, {zj, — 
as depicted in Fig.[T] With separation of pure imaginary 
and true complex values, we denote these N poles as 



-ij s \ s = l,--- ,N r 

-{uj s +i ls )] s = N r + !,-■■ ,N P 



(9) 



together with —z* = lo s — i~/ s in lower plane; thus, N p = 
(N + N r )/2. The parameters 7 S and cu s are all positive. 



The corresponding residues are (Res) 



and (Res) ( _ z .) = e lz ^J{-z* s ) 
identity is inferred from Eq. (JHJ 



= -[(Res) z J*. The last 
Wc have therefore 



C B (t) =^c s e 7st + ^ e 7st [a s cos(uj s t)-b s sm(L0 s t)]. 

(10) 



s=N r + l 



It decomposes Cb(£) into a total of N components. The 
involving coefficients are all real and defined as 

2 4 
c s = — J(-fy a ), a s + ib s = — J[-(oj s + ij s )]. (11) 
ip ip 

The fact that c s is real can be readily proved via its defi- 
nition and Eq. ([5]). Thus, Cg(t) is always a real function. 

The Co(t) component of bath correlation function, 
which is complex in general, arises from the poles of spec- 
tral density J(z). In principle it would not depend on 
PFD/MSD scheme; but in practice it will, for the is- 
sue of consistency to be discussed later. In the following 
HEOM construction, we consider the Drude model, 



J(«) 



rjuj c uj 



(12) 



where r\ is the system-bath coupling strength and uj c the 
cut-off frequency. The Drude spectral density function 
has only one pole in lower plane. It results in 



C {t) = c ( 



(13) 



with 



1 



,-0z 



Z — — 1L0 C 



The above Co(t) expression is exact. However, a consis- 
tency issue arises as the Bose function used in evaluating 
Cb (t) via Eq. (fTTJjl involves the PFD approximation. The 
inconsistency here may be problematic in implementa- 
tion. For example, the divergence of the above exact 
Co(t) at f3us c — 27rm cannot be cancelled out by the ap- 
proximated Cb(£) of Eq. (fTU|) . To overcome this problem, 
we evaluate Co also using PFD of Eq. (J4| with the same 
finite N. 



co 



-irjuj, 



i 



1 



2z 



(14) 



The resulting Co(t) does practically depend on the 
scheme of expansion for Bose function. We shall em- 
phasize the importance of the aforementioned consistent 
treatment between Co(t) and Cb(£)- We have recently 
demonstrated it with different approximation schemes 
[TtJ • We have also confirmed this issue with the present 
PFD scheme, not just for C(t) — Co(t) + Cb(£), but more 
importantly for the numerical performance of the result- 
ing HEOM dynamics to be presented below. 



III. HIERARCHICAL EQUATIONS OF MOTION 
A. HEOM-QDT in PFD scheme 

The generic form of HEOM-QDT reads @ 

p n = -{iC + T n ) Pn + P ir l + pt } + P^- (15) 

Here, C ■ = • ] is the reduced system Liouvillian; 

r n is the decay constant; p n , p„" } , Pn~ } , and p„ +} are well- 
defined auxiliary density operators (ADOs) in the system 
subspace. The reduced system density operator of pri- 
mary interest, defined as p(t) = TrbathPtotai(i), is just 
p{t) ee p Q (t), with pl" } = p^ } = r = 0. The labeling 
index n for a generic ADO p n consists of a set of non- 
negative integers, which are arranged in relation to the 
individual components of bath correlation function in a 
given decomposition scheme. Let n = {n ,n 1; - • • ,njv} 
and no+rti + - • -+nN = n- The latter is used to define the 
tier of p n . The last three terms in Eq. (fT5"|) describes how 
a given p n of n th tier depends on its associated ADOs of 
same tier (plC } ) and neighboring tiers (pn ±} ) , respectively. 

The HEOM QDT formalism is exact and nonpertur- 
bative, assuming only Gaussian bath statistics. It is 
equivalent to the Feynman- Vernon path integral influ- 
ence functional theory of reduced system density oper- 
ator dynamics. Moreover, it also supports the incor- 
poration of the initial system-bath correlation through 
appropriate initial p n ^o(^o) 7^ conditions. HEOM re- 
solves the combined effect of the coupling bath strength 
and memory contents, as they are decomposed, on the 
reduced system dynamics. 

The specific form of theory depends on the way of de- 
composing the bath correlation function C(t). For Drude 
dissipation in the PFD scheme, C(t) [Eq. ([7])] has been 
decomposed into (N + 1) components in Eqs. (fTTJ|) and 
(fT"3f . Therefore, the ADO labeling index n assumes 



n ={?i s =o,i, 



,N. 



}■ (16) 



The composite (iV+1) nonnegative integers are the lead- 
ing orders of individual components of C(t), involved in 
the specified p n . Specifically, no is for the Drude compo- 
nent Co(t) of Eq. (fT3")) . and the other N integers are for 
the N components of Cb(£) of Eq. (p~0|) . respectively. 

In Appendix [Bl we present the standard approach 
to the corresponding HEOM formalism, based on the 
present PFD scheme. The final results are summarized as 
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follows. The parameter T n in Eq. ([15]) collects all damp- 
ing factors. It reads (setting j = lo c ): 

N r N p 

r n = 2J rvy s + ^ (n s + n s )j s . (17) 

s=0 s=JV r +l 

It arises from the derivatives of individual exponential 
components of C (t) , involved in p n . 

The swap term p n " } describes how a given p n depends 
on its associated ADOs of same tier. It reads 



- ^2 ^3 (a>s/b s )\ / n s (n s + l)\b s /a s \ p„ r 

s=N r +l 

- (b s /a s )y / n s (n s + l)\a s /b s \ p nr , (18) 



It arises from the derivatives of the sine and cosine 
components of Cb(£) [Eq. (|10p] involved in p n . The in- 
volving index n^T* differs from n by changing (n s ,n s ) 
to (n s — l,n s + 1) and n^" by changing (n s ,n s ) to 
(n s + l,n s - 1). 

The last two terms in Eq. (fTS"]) describes the tier-down 
and tier-up contributions. They are given by 



Pn — 



~ V n s/\c s \(csQp n7 - c s p n -Q 



s=0 



- Y a s\/n s /\a s \[Q,p n7 ], 

s=N r + l 



(19) 



N r 



n {+} 

Pn 



■^2y/{n s + l)\c 3 \ [Q,p„+] 



s=0 

AT. 



i J2 {V(n s + l)\a s \[Q,Pnt] 
i 

y/(n s + l)\b 3 \ [Q, Pn +]}. 



s=N r + l 



(20) 



They depend explicitly on the system dissipation mode 
Q. The involving index differs from n by changing 
the specified n s to n s ± 1, and Pi f similarly. 

Apparently, the boundary conditions of Pq" } = p^ 1 = 
To = are satisfied. All coefficients/parameters aris- 
ing from the Bose part Cb(£) of bath correlation, 
{a s , 6 S , c s , 7 S , w s } s> o, are real; the last two are positive. 
The Drude damping parameter is set to be 70 = co c . Only 
cq [Eq. fTl|)] is complex. 



ADO if its no + ■ ■ ■ + = h. The number of ADOs at a 
given tier is ^grSr > while the total number of ADOs up 
to level L is 



n=0 



(N + n) \ _ ( L 
Nl n! 



{»}■ 



(21) 



The second identity serves also the definition of {™} 
for abbreviated notion, with the boundary values of 
{ m < } = and {0} = = 1. The working index 
in = 3n -n N = 0, • • • , N - 1 can then be 



N 



n N - 1 



N 



q=0 



}• (22) 



It sorts the index n by tiers and followed by subindices, 



1, 



and 



so that in={o,-,o} = 0) in={i.o,-,o} 

in={0,- ,0,L) =M-\. 

In practice, HEOM should be truncated properly at 
finite level of hierarchy L and decomposition order N . 
By far the truncation of N goes with convergency, but 
that of L are carried out effectively and automatically. 
Apparently, the issue of truncation is closely related to 
central processing unit (CPU) time and memory cost of 
computation. The number Af of total ADOs goes with 
a combinatory law, like the complete configuration inter- 
action treatment in quantum mechanics. Shi et al. have 
recently proposed an efficient, accuracy-controlled, dy- 
namics filtering algorithm [l3|, |3| • The basic idea behind 
is the observation that only a very small fraction of total 
ADOs play roles in HEOM propagation. The filtering 
algorithm sets a specific p n (tj ) = if its matrix elements 
amplitudes are all smaller than the pre-chosen error tol- 
erance. To validate this simple algorithm, all ADOs in- 
cluding the reduced density matrix of the primary inter- 
est should have a uniform filtering-algorithm error toler- 
ance. The present HEOM formalism [Eqs. [fT5 1) -(|2"0 ]) ] has 
been scaled properly to meet this requirement. In con- 
nection with the stochastic field approach to Gaussian- 
Markovian dissipation, the scaled ADOs are just the ex- 
pansion coefficients, with the normalized harmonic wave 
functions used as the basis set for resolving the diffusive 
bath field [T3|. The filtering algorithm keeps only those 
necessary ADOs according to the selected error tolerance. 
Apparently, it also automatically truncates the level of 
hierarchy on-the-fly during numerical propagation. 



IV. NUMERICAL PERFORMANCE AND 
CONCLUDING REMARKS 



B. Remarks on implementation 

To facilitate locating a specified ADO, we also like to 
have a working index scheme to map a set of (TV + 1) 
ordered multiple indices, n = {no,ni, - • • ,wjv}, to an 
integer j n , such that p u = pj n . That p n is an n th -tier 



The performance of the present PFD-based HEOM 
dynamics is exemplified with Fig. for a spin-boson sys- 
tem H = e<j z +Aa x , with the dissipation mode of Q = a z . 
The parameters for the system and coupling bath spec- 
tral density are (e, 77, u> c ) = (1,0.5,5) in unit of A, with 
A = 1660 cm- 1 and T = 298 K. The result of PFD- 
HEOM with TV = 12 in Fig. [3] is practically converged, 
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t 1 i 1 i 1 — r 




t (unit of 2n/A) 

FIG. 3: Time evolution of a spin-boson system density ma- 
trix: (a) the population difference pn(t) — p22(t); (b) the 
deviation of off-diagonal pi2(t) from the exact. The initial 
conditions are all p n (0) = 0; except for p„=o(0) = p(0) where 
pn(0) = 1 and p22(0) = pi2(0) = 0. The specified values 
of (CPU time; A/" e ff) highlight the implementation cost. Espe- 
cially the two underlined values indicate that PFD-HEOM 
converges about 25 times faster than MSD-HEOM. See text 
for details. 

within 0.005 from that with TV = 16. The latter is treated 
as the exact in Fig.[3Jb). The maximum tier of survival 
ADOs is found to be L max — 8, with the filtering algo- 
rithm error tolerance of 10 -7 ; see Sec. lIII Bl or Ref. |13j . 

Included for comparison are also the results of the best 
available MSD-based conventional HEOM theory. This 
reference theory is augmented with the well-established 
Markovian residue correction (MRC) [6j, [7j, in which 
5C(t) = C exac t(i) — Cmsd(£) is approximated by white 
noise for its effect on the MSD-based HEOM dynam- 
ics. Without additional implementation cost, this treat- 
ment significantly reduces the req uired number N for 
converged MSD dynamics @, 0, El EI EE El E3 ■ This 
remarkable MSD-based feature remains true for the sys- 
tem in Fig. [3] However, MRC does not work well with the 
present PFD-HEOM theory. Therefore the comparison 
shown in Fig. [3] is between PFD-HEOM without residue 
correction versus MSD-HEOM with MRC, the best avail- 
able reference as we know. 

Performance is reported in terms of (CPU time; Af c s)- 
The CPU time is by a single Intel(R) Xeon(R) proces- 
sor @3.00 GHz, with the fourth-order Runge-Kutta prop- 
agator and time-step of 0.0015 fs. Another parameter 



Afcs records the largest number of ADOs ever survived 
during propagation with the filtering algorithm. Recall 
that the ADOs in PFD-HEOM is of the maximum sur- 
vival tier level of L max = 8. In contrast, a converged 
MSD-HEOM dynamics requires the level of about 20 
tiers. 

For the purpose of comparison, however, we set L = 
9 in all calculations, including both PFD and MSD 
schemes. Apparently, A4ff in either scheme is orders of 
magnitude smaller than the total number of ADOs, which 
is Af = 497420 for N = 12 or N = 48620 for N = 8, re- 
spectively. The PFD-HEOM is of smaller Af eS , about 
25 ~ 40% of its MSD counterpart with same L and N. 
The enhanced filtering efficiency in PFD scheme here is 
closely related to its complex poles, rather than only pure 
imaginary ones (Fig. [J). The complex poles lead to os- 
cillatory decomposition components of bath correlation 
function [Eq. (fTTJ)) ] , and result in oscillatory cancelation 
in PFD-HEOM dynamics. This right property of the 
PFD scheme may account for the relatively small num- 
ber of survival ADOs after filtering algorithm. 

Performance of PFD-HEOM dynamics should be 
based on the CPU time versus its MSD counterpart of 
the same quality. As mentioned earlier, however, the con- 
verged MSD dynamics for the present system of study is 
too expensive to be worth here. We rather choose for 
assessment by a pair of similar quality but approximated 
results: PFD(7V = 8) versus MSD(7V = 12), with the 
CPU time of 16 minutes versus 397 minutes, respectively, 
see Fig. |31 We have carried out a series of dynamics sim- 
ulations with different system and bath parameters. All 
results show that the performance of PFD-HEOM is su- 
perior over its MSD counterpart by at least an order of 
magnitude. 

In summary, we have constructed the PFD scheme to 
Bose function and HEOM-QDT. The superiority of PFD 
over MSD is apparent. The complex PFD poles lead to 
not just the Bose function expansion more efficient and 
accurate, but also the HEOM construction more com- 
pact. The resulting HEOM dynamics converges with 
smaller L and N, and also accommodates better with 
the propagation filtering algorithm. These factors con- 
tribute the performance of PFD-based HEOM theory, 
which converges a magnitude (about 20 times) faster 
than its MSD counterpart. 
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APPENDIX A: DERIVATION OF BOSE 
FUNCTION IN PFD EXPANSION 

The derivation of Eqs. (j4|) and {5]) starts with the 
Taylor's expansions cosh(y) w J2 n =o 2/ 2n /(2w)! and 
sinh(y) w X) n=0 y 2n+1 /(2n + 1)! for the numerator and 
denominator of the last term in Eq. ((3]), respectively. We 
have then 



where 



B N (y) 



\ + -+B N {x/2). 
2 x 



i Et y 2 "/(2 



2 Eto2/ 2 " +1 /(2«+l)! % 

E^V ra „,2n-l 

n=l (2n+l)! » 



V W 1 7/2" 

Z^n=0 (2n+l)! y 

AT 



IE 



2 " K y 



ii y vt 



(Al) 



(A2) 



By comparing Eqs. (|A5[) and (|A6|) and setting cat = 
a at, we have the following nonzero elements: A^ n — 
-On-i/cn and j4 m)TO+ i = c m /c m+ i for m 7^ TV; others 
are all zeroes. Therefore, 

A m n = <5m+l,n <^miV) (A7) 



with c„ being arbitrary nonzero parameters, except for 
cat = ax- 

In particular, we choose c n = a„_iajv/ajv_i, if every 
a,j =/= 0. The boundary condition of Cat = a at is also 
satisfied. We can recast Eq. (|A7j) as 



(A8) 



It recovers Eq. ([5]) for the polynomial Eq. (|A3|) , where 
a„ = l/(2n+l)!. 



APPENDIX B: DERIVATION OF HEOM-QDT IN 
PFD SCHEME 



with {yj — ±yf£~j\j = 1,- •• , N} being the roots of the 
denominator polynomial; i.e., 



N 

E 



1 



(2n+l) 



■ C = 0. 



(A3) 



The PFD coefficient bj can be determined via 



bf = 2 lim 



= 2 lim 



yB N (yT 



2yE„=i(5^TT)i(yTV^) 



Sn=0 (2n+l)! 

v^N 2n /_ /t7\2n-l 

Z^ n =l (2n+l)! UVWJ 



En=0 (2ri+l)! ("Fy^j) 



(A4) 



The PFD expression of Bose function, Eq. ([J]), is obtained 
by substituting Eqs. ((A2|) and (|A4[) into Eq. (|AT|) . with 
x = 2y. 

To convert the roots of polynomial to eigenvalue prob- 
lem, as Eq. (O, let us consider the general case of 



a + ai£ 



a N £ N = 0. 



(A5) 



We search for the proper matrix A = {A mn }, with the 
eigenvector the form of v = [c%, C2^, ■ • • , cn£ n ~ 1 ] t that 
converts the eigenequation Aw = £v to 



The following derivation of the HEOM in Sec. lIII Al is 
carried out via the well-established Calculus-On-Path- 
Integral Influence-Generating-Functional (COPI-IGF) 
algebra It starts with the path integral influence 
functional for quantum Gaussian dissipation [H, 0, 0] , fol- 
lowed by consecutive time derivatives to resolve in a hier- 
archical manner the involving memory contents 0, & HI • 
Unlike the HEOM theory that can be expressed in op- 
erator level, the path integral formalism has to go with 
representation. Let {\a}} be a generic basis set in the 
system subspace and a = (a, a') for abbreviation, such 
that p(a,t) = p(a,a',t) = (a\p(t)\a'}. Introduce the 
reduced Liouville-space propagator hi via 



p(a,t) = J daoU(a,t;aio,to)p(ac ,to). (Bl) 

Its path-integral expression reads 
/■«[*] 

U(a,t;oi ,t ) = Vae lS[a]/h T[a\e- lS[a]/h . (B2) 

Ja [t ] 

S and T arc the action and influence functionals, respec- 
tively. For Gaussian bath interactions, the latter assumes 
the Feynman- Vernon form that can be recast as [l|, Q 

r[cx}=exp{~^drA[a(r)]B(r;{cx})}, (B3) 



where 



A{a(t)} = Q[a(t)}-Q{a'(t)], 



(B4) 



N 



n-1 



c m £ m ; m = l, 



,N. 



(A6) 



n=l 



B(t; {a}) = -- [B(t; {a}) - B'(t; {a'})] , (B5a) 
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with 



for s = N r + l, 



,N P . 



B(t;{a})= / dTC(t-r)Q[a(T)}, 



(B5b) 



B'(t;{a'})= / dTC*{t- T )Q[a'(T)}. 

J to 



Note that .A[a(t)] depends only on the local time of 
path and assumes readily in operator level as commu- 
tator, A- = [Q, ■}■ In contrast, the functional B(t;{a}) 
contains memory and does not have a simple operator- 
level correspondence. The COPI-IGF algebra provides a 
proper hierarchy to resolve the memory contents involved 
in B(t;{a}). 

To proceed, we decompose Bit; {a}) according to the 
decomposition of bath correlation function C (t) [Eq. ([7]) 
with Eqs. ^TDJ) and (US])]. We have 



N r 



(B6) 



s=0 s=N r +l 

where (noting that cq is complex while others are real) 



B s = --(c s B s -c* s B' s ); a = 0,l,-- ,N r , (B7a) 
B s = ~a s (B s - B' s ); s = N r + 1, • • • , N p , (B7b) 
B s = l -b s {B s - B' s ); s = N r + 1, ■ ■ ■ , N p . (B7c) 



with (denoting uj s <N r = 0) 



B s = / dr e -7.(t-^) CO s[uj s (t - T)]Q[a(r)], 



B s = / dre-^(*- T )sin[ Ws (i-r)]Q[a(r)]. 



(B8) 



The time derivative of Eq. (|B8|1 reads 

d t B s = - ls B s - uj s B s + Q[a(t)}, 



(B9) 



d t B s = —%B S + uj s B s . 
Thus the memory functionals of Eq. (|B7[) satisfy 

d t B s = - ls B s - l -{c s Q[a{t)] - c* s Q[a'(t)]}, (BlOa) 

for s = 0, 1, • • • ,N r ; while 

d t B s = - ls B s + ^u s B s - % -a s A[ot{t)}, (BlOb) 



d t B s = ~ j s B s uj s B s , 

a., 



(BlOc) 



The above equations of motion for {B s , B s } are closed, 
with the inhomogeneous terms depending only on the lo- 
cal time. This is right the property for these {B S ,B S } 
defined in Eq. (|B7[) to be the IGFs for the desired hierar- 
chy construction. The auxiliary influence functionals T n , 
with the labeling index n of Eq. (116|) . are now obtained 
via the IGFs as [9| 



/cr„ 



s=N r + l 



(Bll) 



which specifies the ADO as p„(t) = U n (t, to)p(to) with 



/■«[*] 

U n (a,t;a ,to) = / Dae iSW/s f n [a] e -' SH/t 

J a [t ] 



(B12) 



Included in Eq. (|B11[) is also a proper scaling param- 
eter o n . This is for the pur pos e of applying simple and 
efficient filtering algorithm [13( in HEOM propagation. 
This parameter scales the specified p n to be of not just 
the same unit, but also about the same error tolerance as 
the reduced system densit y op erator p = po of primary 
interest. It is given by [l3l.ll4|: 



<r n = l[(\c s \ n °n s \) ■ J] (\a s \ n 'n s \\b s \ n >nj). (B13) 



s=Q 



s=N r + l 



The HEOM [Eqs. (JIH])— (EOD] can now be obtained read- 
ily by taking the time derivative on p n (t) or its propaga- 
tor U n of Eq. (|B12p . Various terms in Eqs. |T5]) and (fT7| - 
(I20p result from the derivatives on various components in 
Eqs. (|Bllj) and Eq. (|B12|) . as detailed as follows. The co- 
herent Liouvillian dynamics part in Eq. (|15p arises from 
the derivative of the action functionals in Eq. (|B12p ; The 
tier-up p { rt } term [Eq. (20])] is from d t T = -^A{BT), fol- 
lowing by using Eqs. (|B6p and (jBlip . and the operator- 
level form of A- = [Q, •]; Finally, T n , p„ , and p„~ } 
[Eqs. (fT7|) - (fT9)) ] collect the decay, swap, and inhomoge- 
neous terms of Eqs. (|B10|) . respectively. Note that the 
scaling parameters, tr n of Eq. (|B13|) . are also involved. 
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